## Section 3.2.3 of the pre-analysis plan: "Overall Robustness to Data
## Collection Procedures: Sensitivity Normalized by Precision"
##
## Kevin Quinn
## University of Michigan
##
## 8/7/2019
##

library(TopKLists)
library(coda)

set.seed(77810)

sink(file="MM3.2.3.output.txt")



########################################################################
## START a1-b1 comparisons
########################################################################

## load MCMC output
load("MM3.2.a1.M1.out.Rda")
a1.M1.out <- a1.M1.out[sample(1:nrow(a1.M1.out), nrow(a1.M1.out),
                              replace=FALSE), ]
load("MM3.2.a1.M2.out.Rda")
a1.M2.out <- a1.M2.out[sample(1:nrow(a1.M2.out), nrow(a1.M2.out),
                              replace=FALSE), ]
load("MM3.2.a1.pairwise.out.Rda")
a1.pairwise.out <- a1.pairwise.out[sample(1:nrow(a1.pairwise.out),
                                          nrow(a1.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b1.M1.out.Rda")
b1.M1.out <- b1.M1.out[sample(1:nrow(b1.M1.out), nrow(b1.M1.out),
                              replace=FALSE), ]

load("MM3.2.b1.M2.out.Rda")
b1.M2.out <- b1.M2.out[sample(1:nrow(b1.M2.out), nrow(b1.M2.out),
                              replace=FALSE), ]
load("MM3.2.b1.pairwise.out.Rda")
b1.pairwise.out <- b1.pairwise.out[sample(1:nrow(b1.pairwise.out),
                                          nrow(b1.pairwise.out),
                              replace=FALSE), ]


## make MCMC output comparable
a1.M1.out <- a1.M1.out[,1:48]
a1.M2.out <- a1.M2.out[,1:48]
b1.M1.out <- b1.M1.out[,1:48]
b1.M2.out <- b1.M2.out[,1:48]

colnames(a1.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a1.M1.out))
colnames(a1.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a1.M2.out))
colnames(b1.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b1.M1.out))
colnames(b1.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b1.M2.out))
colnames(a1.pairwise.out) <- gsub("theta\\.", "", colnames(a1.pairwise.out))
colnames(b1.pairwise.out) <- gsub("theta\\.", "", colnames(b1.pairwise.out))

for (j in 2:48){
    a1.M1.out[,j] <- a1.M1.out[,1] + a1.M1.out[,j]
    a1.M2.out[,j] <- a1.M2.out[,1] + a1.M2.out[,j]
    b1.M1.out[,j] <- b1.M1.out[,1] + b1.M1.out[,j]
    b1.M2.out[,j] <- b1.M2.out[,1] + b1.M2.out[,j]
}
colnames(a1.M1.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(a1.M2.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(b1.M1.out)[1] <- colnames(a1.pairwise.out)[1]
colnames(b1.M2.out)[1] <- colnames(a1.pairwise.out)[1]


load("MM3.2.2.D.hat.a1-b1.data.Rda")
denom.spear.M1 <- (mean(D.hat.data$D.spear.a1.M1) +
                   mean(D.hat.data$D.spear.b1.M1))/2
denom.spear.M2 <- (mean(D.hat.data$D.spear.a1.M2) +
                   mean(D.hat.data$D.spear.b1.M2))/2
denom.spear.pair <- (mean(D.hat.data$D.spear.a1.pair) +
                     mean(D.hat.data$D.spear.b1.pair))/2

denom.kend.M1 <- (mean(D.hat.data$D.kend.a1.M1) +
                   mean(D.hat.data$D.kend.b1.M1))/2
denom.kend.M2 <- (mean(D.hat.data$D.kend.a1.M2) +
                   mean(D.hat.data$D.kend.b1.M2))/2
denom.kend.pair <- (mean(D.hat.data$D.kend.a1.pair) +
                     mean(D.hat.data$D.kend.b1.pair))/2


M <- nrow(a1.M1.out)


## normalized distance between ranks
d.spear.M1 <- rep(NA, M)
d.spear.M2 <- rep(NA, M)
d.spear.pair <- rep(NA, M)

d.kend.M1 <- rep(NA, M)
d.kend.M2 <- rep(NA, M)
d.kend.pair <- rep(NA, M)

spear.M1.pair <- rep(NA, M)
kend.M1.pair <- rep(NA, M)
spear.M2.pair <- rep(NA, M)
kend.M2.pair <- rep(NA, M)
for (iter in 1:M){
    d.spear.M1[iter] <- Spearman(rank(a1.M1.out[iter,]),
                                 rank(b1.M1.out[iter,]),
                                 48, 48) / denom.spear.M1
    d.spear.M2[iter] <- Spearman(rank(a1.M2.out[iter,]),
                                 rank(b1.M2.out[iter,]),
                                 48, 48) / denom.spear.M2
    d.spear.pair[iter] <- Spearman(rank(a1.pairwise.out[iter,]),
                                   rank(b1.pairwise.out[iter,]),
                                   48, 48) / denom.spear.pair
    
    d.kend.M1[iter] <- Kendall2Lists(rank(a1.M1.out[iter,]),
                                     rank(b1.M1.out[iter,]),
                                 48, 48, 48) / denom.kend.M1
    d.kend.M2[iter] <- Kendall2Lists(rank(a1.M2.out[iter,]),
                                     rank(b1.M2.out[iter,]),
                                 48, 48, 48) / denom.kend.M2
    d.kend.pair[iter] <- Kendall2Lists(rank(a1.pairwise.out[iter,]),
                                       rank(b1.pairwise.out[iter,]),
                                 48, 48, 48) / denom.kend.pair


    spear.M1.pair[iter] <- d.spear.pair[iter] < d.spear.M1[iter]
    spear.M2.pair[iter] <- d.spear.pair[iter] < d.spear.M2[iter]
    kend.M1.pair[iter] <- d.kend.pair[iter] < d.kend.M1[iter]
    kend.M2.pair[iter] <- d.kend.pair[iter] < d.kend.M2[iter]
}




cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat(" Normalized distances between ranks (a1 vs. b1)\n")
cat("-----------------------------------------------------------\n")
cat("d.spear.M1 =  ", mean(d.spear.M1), "    [", quantile(d.spear.M1, probs=0.025), ", ", quantile(d.spear.M1, probs=0.975), "]\n", sep="")
cat("d.spear.M2 =  ", mean(d.spear.M2), "    [", quantile(d.spear.M2, probs=0.025), ", ", quantile(d.spear.M2, probs=0.975), "]\n", sep="")
cat("d.spear.pair =", mean(d.spear.pair), "    [", quantile(d.spear.pair, probs=0.025), ", ", quantile(d.spear.pair, probs=0.975), "]\n", sep="")
cat("\n")
cat("d.kend.M1 =  ", mean(d.kend.M1), "    [", quantile(d.kend.M1, probs=0.025), ", ", quantile(d.kend.M1, probs=0.975), "]\n", sep="")
cat("d.kend.M2 =  ", mean(d.kend.M2), "    [", quantile(d.kend.M2, probs=0.025), ", ", quantile(d.kend.M2, probs=0.975), "]\n", sep="")
cat("d.kend.pair =", mean(d.kend.pair), "    [", quantile(d.kend.pair, probs=0.025), ", ", quantile(d.kend.pair, probs=0.975), "]\n", sep="")

cat("===========================================================\n")
cat("a1 - b1 normalized comparisons (coder race corr. with photos)\n")
cat("-----------------------------------------------------------\n")
cat("M1 vs. pairwise: \n")
cat("  Spearman's Footrule\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(spear.M1.pair), "\n")
cat("  Kendall's Tau\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(kend.M1.pair), "\n")
cat("-----------------------------------------------------------\n")
cat("M2 vs. pairwise: \n")
cat("  Spearman's Footrule\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(spear.M2.pair), "\n")
cat("  Kendall's Tau\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(kend.M2.pair), "\n")
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")






########################################################################
## END a1-b1 comparisons
########################################################################












########################################################################
## START a2-b2 comparisons
########################################################################

## load MCMC output
load("MM3.2.a2.M1.out.Rda")
a2.M1.out <- a2.M1.out[sample(1:nrow(a2.M1.out), nrow(a2.M1.out),
                              replace=FALSE), ]
load("MM3.2.a2.M2.out.Rda")
a2.M2.out <- a2.M2.out[sample(1:nrow(a2.M2.out), nrow(a2.M2.out),
                              replace=FALSE), ]
load("MM3.2.a2.pairwise.out.Rda")
a2.pairwise.out <- a2.pairwise.out[sample(1:nrow(a2.pairwise.out),
                                          nrow(a2.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b2.M1.out.Rda")
b2.M1.out <- b2.M1.out[sample(1:nrow(b2.M1.out), nrow(b2.M1.out),
                              replace=FALSE), ]

load("MM3.2.b2.M2.out.Rda")
b2.M2.out <- b2.M2.out[sample(1:nrow(b2.M2.out), nrow(b2.M2.out),
                              replace=FALSE), ]
load("MM3.2.b2.pairwise.out.Rda")
b2.pairwise.out <- b2.pairwise.out[sample(1:nrow(b2.pairwise.out),
                                          nrow(b2.pairwise.out),
                              replace=FALSE), ]


## make MCMC output comparable
a2.M1.out <- a2.M1.out[,1:48]
a2.M2.out <- a2.M2.out[,1:48]
b2.M1.out <- b2.M1.out[,1:48]
b2.M2.out <- b2.M2.out[,1:48]

colnames(a2.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a2.M1.out))
colnames(a2.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a2.M2.out))
colnames(b2.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b2.M1.out))
colnames(b2.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b2.M2.out))
colnames(a2.pairwise.out) <- gsub("theta\\.", "", colnames(a2.pairwise.out))
colnames(b2.pairwise.out) <- gsub("theta\\.", "", colnames(b2.pairwise.out))

for (j in 2:48){
    a2.M1.out[,j] <- a2.M1.out[,1] + a2.M1.out[,j]
    a2.M2.out[,j] <- a2.M2.out[,1] + a2.M2.out[,j]
    b2.M1.out[,j] <- b2.M1.out[,1] + b2.M1.out[,j]
    b2.M2.out[,j] <- b2.M2.out[,1] + b2.M2.out[,j]
}
colnames(a2.M1.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(a2.M2.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(b2.M1.out)[1] <- colnames(a2.pairwise.out)[1]
colnames(b2.M2.out)[1] <- colnames(a2.pairwise.out)[1]


load("MM3.2.2.D.hat.a2-b2.data.Rda")
denom.spear.M1 <- (mean(D.hat.data$D.spear.a2.M1) +
                   mean(D.hat.data$D.spear.b2.M1))/2
denom.spear.M2 <- (mean(D.hat.data$D.spear.a2.M2) +
                   mean(D.hat.data$D.spear.b2.M2))/2
denom.spear.pair <- (mean(D.hat.data$D.spear.a2.pair) +
                     mean(D.hat.data$D.spear.b2.pair))/2

denom.kend.M1 <- (mean(D.hat.data$D.kend.a2.M1) +
                   mean(D.hat.data$D.kend.b2.M1))/2
denom.kend.M2 <- (mean(D.hat.data$D.kend.a2.M2) +
                   mean(D.hat.data$D.kend.b2.M2))/2
denom.kend.pair <- (mean(D.hat.data$D.kend.a2.pair) +
                     mean(D.hat.data$D.kend.b2.pair))/2


M <- nrow(a2.M1.out)


## normalized distance between ranks
d.spear.M1 <- rep(NA, M)
d.spear.M2 <- rep(NA, M)
d.spear.pair <- rep(NA, M)

d.kend.M1 <- rep(NA, M)
d.kend.M2 <- rep(NA, M)
d.kend.pair <- rep(NA, M)

spear.M1.pair <- rep(NA, M)
kend.M1.pair <- rep(NA, M)
spear.M2.pair <- rep(NA, M)
kend.M2.pair <- rep(NA, M)
for (iter in 1:M){
    d.spear.M1[iter] <- Spearman(rank(a2.M1.out[iter,]),
                                 rank(b2.M1.out[iter,]),
                                 48, 48) / denom.spear.M1
    d.spear.M2[iter] <- Spearman(rank(a2.M2.out[iter,]),
                                 rank(b2.M2.out[iter,]),
                                 48, 48) / denom.spear.M2
    d.spear.pair[iter] <- Spearman(rank(a2.pairwise.out[iter,]),
                                   rank(b2.pairwise.out[iter,]),
                                   48, 48) / denom.spear.pair
    
    d.kend.M1[iter] <- Kendall2Lists(rank(a2.M1.out[iter,]),
                                 rank(b2.M1.out[iter,]),
                                 48, 48, 48) / denom.kend.M1
    d.kend.M2[iter] <- Kendall2Lists(rank(a2.M2.out[iter,]),
                                 rank(b2.M2.out[iter,]),
                                 48, 48, 48) / denom.kend.M2
    d.kend.pair[iter] <- Kendall2Lists(rank(a2.pairwise.out[iter,]),
                                 rank(b2.pairwise.out[iter,]),
                                 48, 48, 48) / denom.kend.pair


    spear.M1.pair[iter] <- d.spear.pair[iter] < d.spear.M1[iter]
    spear.M2.pair[iter] <- d.spear.pair[iter] < d.spear.M2[iter]
    kend.M1.pair[iter] <- d.kend.pair[iter] < d.kend.M1[iter]
    kend.M2.pair[iter] <- d.kend.pair[iter] < d.kend.M2[iter]
}




cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat(" Normalized distances between ranks (a2 vs. b2)\n")
cat("-----------------------------------------------------------\n")
cat("d.spear.M1 =  ", mean(d.spear.M1), "    [", quantile(d.spear.M1, probs=0.025), ", ", quantile(d.spear.M1, probs=0.975), "]\n", sep="")
cat("d.spear.M2 =  ", mean(d.spear.M2), "    [", quantile(d.spear.M2, probs=0.025), ", ", quantile(d.spear.M2, probs=0.975), "]\n", sep="")
cat("d.spear.pair =", mean(d.spear.pair), "    [", quantile(d.spear.pair, probs=0.025), ", ", quantile(d.spear.pair, probs=0.975), "]\n", sep="")
cat("\n")
cat("d.kend.M1 =  ", mean(d.kend.M1), "    [", quantile(d.kend.M1, probs=0.025), ", ", quantile(d.kend.M1, probs=0.975), "]\n", sep="")
cat("d.kend.M2 =  ", mean(d.kend.M2), "    [", quantile(d.kend.M2, probs=0.025), ", ", quantile(d.kend.M2, probs=0.975), "]\n", sep="")
cat("d.kend.pair =", mean(d.kend.pair), "    [", quantile(d.kend.pair, probs=0.025), ", ", quantile(d.kend.pair, probs=0.975), "]\n", sep="")

cat("===========================================================\n")
cat("a2 - b2 normalized comparisons (distractor race corr w. photos)\n")
cat("-----------------------------------------------------------\n")
cat("M1 vs. pairwise: \n")
cat("  Spearman's Footrule\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(spear.M1.pair), "\n")
cat("  Kendall's Tau\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(kend.M1.pair), "\n")
cat("-----------------------------------------------------------\n")
cat("M2 vs. pairwise: \n")
cat("  Spearman's Footrule\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(spear.M2.pair), "\n")
cat("  Kendall's Tau\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(kend.M2.pair), "\n")
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")






########################################################################
## END a2-b2 comparisons
########################################################################
















########################################################################
## START a3-b3 comparisons
########################################################################

## load MCMC output
load("MM3.2.a3.M1.out.Rda")
a3.M1.out <- a3.M1.out[sample(1:nrow(a3.M1.out), nrow(a3.M1.out),
                              replace=FALSE), ]
load("MM3.2.a3.M2.out.Rda")
a3.M2.out <- a3.M2.out[sample(1:nrow(a3.M2.out), nrow(a3.M2.out),
                              replace=FALSE), ]
load("MM3.2.a3.pairwise.out.Rda")
a3.pairwise.out <- a3.pairwise.out[sample(1:nrow(a3.pairwise.out),
                                          nrow(a3.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b3.M1.out.Rda")
b3.M1.out <- b3.M1.out[sample(1:nrow(b3.M1.out), nrow(b3.M1.out),
                              replace=FALSE), ]

load("MM3.2.b3.M2.out.Rda")
b3.M2.out <- b3.M2.out[sample(1:nrow(b3.M2.out), nrow(b3.M2.out),
                              replace=FALSE), ]
load("MM3.2.b3.pairwise.out.Rda")
b3.pairwise.out <- b3.pairwise.out[sample(1:nrow(b3.pairwise.out),
                                          nrow(b3.pairwise.out),
                              replace=FALSE), ]


## make MCMC output comparable
a3.M1.out <- a3.M1.out[,1:48]
a3.M2.out <- a3.M2.out[,1:48]
b3.M1.out <- b3.M1.out[,1:48]
b3.M2.out <- b3.M2.out[,1:48]

colnames(a3.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a3.M1.out))
colnames(a3.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a3.M2.out))
colnames(b3.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b3.M1.out))
colnames(b3.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b3.M2.out))
colnames(a3.pairwise.out) <- gsub("theta\\.", "", colnames(a3.pairwise.out))
colnames(b3.pairwise.out) <- gsub("theta\\.", "", colnames(b3.pairwise.out))

for (j in 2:48){
    a3.M1.out[,j] <- a3.M1.out[,1] + a3.M1.out[,j]
    a3.M2.out[,j] <- a3.M2.out[,1] + a3.M2.out[,j]
    b3.M1.out[,j] <- b3.M1.out[,1] + b3.M1.out[,j]
    b3.M2.out[,j] <- b3.M2.out[,1] + b3.M2.out[,j]
}
colnames(a3.M1.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(a3.M2.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(b3.M1.out)[1] <- colnames(a3.pairwise.out)[1]
colnames(b3.M2.out)[1] <- colnames(a3.pairwise.out)[1]


load("MM3.2.2.D.hat.a3-b3.data.Rda")
denom.spear.M1 <- (mean(D.hat.data$D.spear.a3.M1) +
                   mean(D.hat.data$D.spear.b3.M1))/2
denom.spear.M2 <- (mean(D.hat.data$D.spear.a3.M2) +
                   mean(D.hat.data$D.spear.b3.M2))/2
denom.spear.pair <- (mean(D.hat.data$D.spear.a3.pair) +
                     mean(D.hat.data$D.spear.b3.pair))/2

denom.kend.M1 <- (mean(D.hat.data$D.kend.a3.M1) +
                   mean(D.hat.data$D.kend.b3.M1))/2
denom.kend.M2 <- (mean(D.hat.data$D.kend.a3.M2) +
                   mean(D.hat.data$D.kend.b3.M2))/2
denom.kend.pair <- (mean(D.hat.data$D.kend.a3.pair) +
                     mean(D.hat.data$D.kend.b3.pair))/2


M <- nrow(a3.M1.out)


## normalized distance between ranks
d.spear.M1 <- rep(NA, M)
d.spear.M2 <- rep(NA, M)
d.spear.pair <- rep(NA, M)

d.kend.M1 <- rep(NA, M)
d.kend.M2 <- rep(NA, M)
d.kend.pair <- rep(NA, M)

spear.M1.pair <- rep(NA, M)
kend.M1.pair <- rep(NA, M)
spear.M2.pair <- rep(NA, M)
kend.M2.pair <- rep(NA, M)
for (iter in 1:M){
    d.spear.M1[iter] <- Spearman(rank(a3.M1.out[iter,]),
                                 rank(b3.M1.out[iter,]),
                                 48, 48) / denom.spear.M1
    d.spear.M2[iter] <- Spearman(rank(a3.M2.out[iter,]),
                                 rank(b3.M2.out[iter,]),
                                 48, 48) / denom.spear.M2
    d.spear.pair[iter] <- Spearman(rank(a3.pairwise.out[iter,]),
                                   rank(b3.pairwise.out[iter,]),
                                   48, 48) / denom.spear.pair
    
    d.kend.M1[iter] <- Kendall2Lists(rank(a3.M1.out[iter,]),
                                 rank(b3.M1.out[iter,]),
                                 48, 48, 48) / denom.kend.M1
    d.kend.M2[iter] <- Kendall2Lists(rank(a3.M2.out[iter,]),
                                 rank(b3.M2.out[iter,]),
                                 48, 48, 48) / denom.kend.M2
    d.kend.pair[iter] <- Kendall2Lists(rank(a3.pairwise.out[iter,]),
                                 rank(b3.pairwise.out[iter,]),
                                 48, 48, 48) / denom.kend.pair


    spear.M1.pair[iter] <- d.spear.pair[iter] < d.spear.M1[iter]
    spear.M2.pair[iter] <- d.spear.pair[iter] < d.spear.M2[iter]
    kend.M1.pair[iter] <- d.kend.pair[iter] < d.kend.M1[iter]
    kend.M2.pair[iter] <- d.kend.pair[iter] < d.kend.M2[iter]
}




cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat(" Normalized distances between ranks (a3 vs. b3)\n")
cat("-----------------------------------------------------------\n")
cat("d.spear.M1 =  ", mean(d.spear.M1), "    [", quantile(d.spear.M1, probs=0.025), ", ", quantile(d.spear.M1, probs=0.975), "]\n", sep="")
cat("d.spear.M2 =  ", mean(d.spear.M2), "    [", quantile(d.spear.M2, probs=0.025), ", ", quantile(d.spear.M2, probs=0.975), "]\n", sep="")
cat("d.spear.pair =", mean(d.spear.pair), "    [", quantile(d.spear.pair, probs=0.025), ", ", quantile(d.spear.pair, probs=0.975), "]\n", sep="")
cat("\n")
cat("d.kend.M1 =  ", mean(d.kend.M1), "    [", quantile(d.kend.M1, probs=0.025), ", ", quantile(d.kend.M1, probs=0.975), "]\n", sep="")
cat("d.kend.M2 =  ", mean(d.kend.M2), "    [", quantile(d.kend.M2, probs=0.025), ", ", quantile(d.kend.M2, probs=0.975), "]\n", sep="")
cat("d.kend.pair =", mean(d.kend.pair), "    [", quantile(d.kend.pair, probs=0.025), ", ", quantile(d.kend.pair, probs=0.975), "]\n", sep="")

cat("===========================================================\n")
cat("a3 - b3 normalized comparisons (black vs. white coders)\n")
cat("-----------------------------------------------------------\n")
cat("M1 vs. pairwise: \n")
cat("  Spearman's Footrule\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(spear.M1.pair), "\n")
cat("  Kendall's Tau\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(kend.M1.pair), "\n")
cat("-----------------------------------------------------------\n")
cat("M2 vs. pairwise: \n")
cat("  Spearman's Footrule\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(spear.M2.pair), "\n")
cat("  Kendall's Tau\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(kend.M2.pair), "\n")
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")






########################################################################
## END a3-b3 comparisons
########################################################################










########################################################################
## START a4-b4 comparisons
########################################################################

## load MCMC output
load("MM3.2.a4.M1.out.Rda")
a4.M1.out <- a4.M1.out[sample(1:nrow(a4.M1.out), nrow(a4.M1.out),
                              replace=FALSE), ]
load("MM3.2.a4.M2.out.Rda")
a4.M2.out <- a4.M2.out[sample(1:nrow(a4.M2.out), nrow(a4.M2.out),
                              replace=FALSE), ]
load("MM3.2.a4.pairwise.out.Rda")
a4.pairwise.out <- a4.pairwise.out[sample(1:nrow(a4.pairwise.out),
                                          nrow(a4.pairwise.out),
                              replace=FALSE), ]

load("MM3.2.b4.M1.out.Rda")
b4.M1.out <- b4.M1.out[sample(1:nrow(b4.M1.out), nrow(b4.M1.out),
                              replace=FALSE), ]

load("MM3.2.b4.M2.out.Rda")
b4.M2.out <- b4.M2.out[sample(1:nrow(b4.M2.out), nrow(b4.M2.out),
                              replace=FALSE), ]
load("MM3.2.b4.pairwise.out.Rda")
b4.pairwise.out <- b4.pairwise.out[sample(1:nrow(b4.pairwise.out),
                                          nrow(b4.pairwise.out),
                              replace=FALSE), ]


## make MCMC output comparable
a4.M1.out <- a4.M1.out[,1:48]
a4.M2.out <- a4.M2.out[,1:48]
b4.M1.out <- b4.M1.out[,1:48]
b4.M2.out <- b4.M2.out[,1:48]

colnames(a4.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a4.M1.out))
colnames(a4.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(a4.M2.out))
colnames(b4.M1.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b4.M1.out))
colnames(b4.M2.out) <- gsub("as.factor\\(photoID\\)", "", colnames(b4.M2.out))
colnames(a4.pairwise.out) <- gsub("theta\\.", "", colnames(a4.pairwise.out))
colnames(b4.pairwise.out) <- gsub("theta\\.", "", colnames(b4.pairwise.out))

for (j in 2:48){
    a4.M1.out[,j] <- a4.M1.out[,1] + a4.M1.out[,j]
    a4.M2.out[,j] <- a4.M2.out[,1] + a4.M2.out[,j]
    b4.M1.out[,j] <- b4.M1.out[,1] + b4.M1.out[,j]
    b4.M2.out[,j] <- b4.M2.out[,1] + b4.M2.out[,j]
}
colnames(a4.M1.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(a4.M2.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(b4.M1.out)[1] <- colnames(a4.pairwise.out)[1]
colnames(b4.M2.out)[1] <- colnames(a4.pairwise.out)[1]


load("MM3.2.2.D.hat.a4-b4.data.Rda")
denom.spear.M1 <- (mean(D.hat.data$D.spear.a4.M1) +
                   mean(D.hat.data$D.spear.b4.M1))/2
denom.spear.M2 <- (mean(D.hat.data$D.spear.a4.M2) +
                   mean(D.hat.data$D.spear.b4.M2))/2
denom.spear.pair <- (mean(D.hat.data$D.spear.a4.pair) +
                     mean(D.hat.data$D.spear.b4.pair))/2

denom.kend.M1 <- (mean(D.hat.data$D.kend.a4.M1) +
                   mean(D.hat.data$D.kend.b4.M1))/2
denom.kend.M2 <- (mean(D.hat.data$D.kend.a4.M2) +
                   mean(D.hat.data$D.kend.b4.M2))/2
denom.kend.pair <- (mean(D.hat.data$D.kend.a4.pair) +
                     mean(D.hat.data$D.kend.b4.pair))/2


M <- nrow(a4.M1.out)


## normalized distance between ranks
d.spear.M1 <- rep(NA, M)
d.spear.M2 <- rep(NA, M)
d.spear.pair <- rep(NA, M)

d.kend.M1 <- rep(NA, M)
d.kend.M2 <- rep(NA, M)
d.kend.pair <- rep(NA, M)

spear.M1.pair <- rep(NA, M)
kend.M1.pair <- rep(NA, M)
spear.M2.pair <- rep(NA, M)
kend.M2.pair <- rep(NA, M)
for (iter in 1:M){
    d.spear.M1[iter] <- Spearman(rank(a4.M1.out[iter,]),
                                 rank(b4.M1.out[iter,]),
                                 48, 48) / denom.spear.M1
    d.spear.M2[iter] <- Spearman(rank(a4.M2.out[iter,]),
                                 rank(b4.M2.out[iter,]),
                                 48, 48) / denom.spear.M2
    d.spear.pair[iter] <- Spearman(rank(a4.pairwise.out[iter,]),
                                   rank(b4.pairwise.out[iter,]),
                                   48, 48) / denom.spear.pair
    
    d.kend.M1[iter] <- Kendall2Lists(rank(a4.M1.out[iter,]),
                                 rank(b4.M1.out[iter,]),
                                 48, 48, 48) / denom.kend.M1
    d.kend.M2[iter] <- Kendall2Lists(rank(a4.M2.out[iter,]),
                                 rank(b4.M2.out[iter,]),
                                 48, 48, 48) / denom.kend.M2
    d.kend.pair[iter] <- Kendall2Lists(rank(a4.pairwise.out[iter,]),
                                 rank(b4.pairwise.out[iter,]),
                                 48, 48, 48) / denom.kend.pair


    spear.M1.pair[iter] <- d.spear.pair[iter] < d.spear.M1[iter]
    spear.M2.pair[iter] <- d.spear.pair[iter] < d.spear.M2[iter]
    kend.M1.pair[iter] <- d.kend.pair[iter] < d.kend.M1[iter]
    kend.M2.pair[iter] <- d.kend.pair[iter] < d.kend.M2[iter]
}




cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat(" Normalized distances between ranks (a4 vs. b4)\n")
cat("-----------------------------------------------------------\n")
cat("d.spear.M1 =  ", mean(d.spear.M1), "    [", quantile(d.spear.M1, probs=0.025), ", ", quantile(d.spear.M1, probs=0.975), "]\n", sep="")
cat("d.spear.M2 =  ", mean(d.spear.M2), "    [", quantile(d.spear.M2, probs=0.025), ", ", quantile(d.spear.M2, probs=0.975), "]\n", sep="")
cat("d.spear.pair =", mean(d.spear.pair), "    [", quantile(d.spear.pair, probs=0.025), ", ", quantile(d.spear.pair, probs=0.975), "]\n", sep="")
cat("\n")
cat("d.kend.M1 =  ", mean(d.kend.M1), "    [", quantile(d.kend.M1, probs=0.025), ", ", quantile(d.kend.M1, probs=0.975), "]\n", sep="")
cat("d.kend.M2 =  ", mean(d.kend.M2), "    [", quantile(d.kend.M2, probs=0.025), ", ", quantile(d.kend.M2, probs=0.975), "]\n", sep="")
cat("d.kend.pair =", mean(d.kend.pair), "    [", quantile(d.kend.pair, probs=0.025), ", ", quantile(d.kend.pair, probs=0.975), "]\n", sep="")

cat("===========================================================\n")
cat("a4 - b4 normalized comparisons (black vs. white distractors)\n")
cat("-----------------------------------------------------------\n")
cat("M1 vs. pairwise: \n")
cat("  Spearman's Footrule\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(spear.M1.pair), "\n")
cat("  Kendall's Tau\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(kend.M1.pair), "\n")
cat("-----------------------------------------------------------\n")
cat("M2 vs. pairwise: \n")
cat("  Spearman's Footrule\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(spear.M2.pair), "\n")
cat("  Kendall's Tau\n")
cat("    Pr(d(theta^a, theta^b)/D < d(beta^a, beta^b)/D) =", mean(kend.M2.pair), "\n")
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n")






########################################################################
## END a4-b4 comparisons
########################################################################


sink()
